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Abstract 

Wc perform a numerical study of black hole formation from the spherically symmetric col- 
lapse of a massless scalar field. The calculations are done in Painleve-GuUstrand (PG) coordi- 
nates that extend across apparent horizons and allow the numerical evolution to proceed until 
the onset of singularity formation. We generate spacetime maps of the collapse and illustrate the 
evolution of apparent horizons and trapping surfaces for various initial data. We also study the 
critical behaviour and find the expected Choptuik scaling with universal values for the critical 
exponent and echoing period consistent with the accepted values of 7 = 0.374 and A = 3.44, 
respectively. The subcritical curvature scaling exhibits the expected oscillatory behaviour but 
the form of the periodic oscillations in the supercritical mass scaling relation, while universal 
with respect to initial PG data, is non-standard: it is non-sinusoidal with large amplitude cusps. 



1 Introduction 



Black holes are among the most interesting solutions to Einstein's field equations. They demonstrate 
in various ways the incompleteness of our classical understanding of gravity and provide hints 
toward developing a quantum theory. The relevant issues are the endpoint of Hawking radiation, 
the breakdown of general relativity at curvature singularities and the microscopic source of black 
hole entropy. Given the lack of experimental guidance, black holes provide a vital testing ground 
for new ideas in quantum gravity. 

Despite a great deal of analytic and numerical work, relatively little is known about the dynamics 
of black hole formation except in special, highly symmetric situations, since the equations for 
collapsing matter fields tend to be very difficult to work with. Important pioneering work in this 
area was done by Choptuik [1] in the early 1990's, who took advantage of modern computational 
methods to numerically create small mass black holes from the collapse of spherically symmetric 
massless scalar fields. In the small mass limit he observed unexpected critical behaviour that has 
since been confirmed in numerous papers for a variety of matter fields (see |2j for an in-depth 
review) . 

There are three fundamental properties of the critical collapse that are relevant here. First, 
there is a critical zero mass black hole solution that plays the role of an intermediate attractor for 
the dynamical evolution. In the case of massless scalar fields the critical solution exhibits discrete 
self similarity (DSS) described in terms of a generic set of field variables by 

Z,(r,t,) = Z,(e"^^r,e'"^i,), (1.1) 

where is the critical solution, r and tg are the Schwarzshild radius and time, m is some integer 
> and A is a dimensionless number. To explain this more clearly, consider the profiles of the 
critical solution field variables frozen in time at tsg. After a certain amount of time 6ts elapses, the 
same profiles will be found on a scale that is smaller by a factor of e^. After an additional time 
6ts/e^ has passed, the same profile will again be found on a scale that is e^^ times smaller than at 
tsQ. This gives rise to the more descriptive name for DSS in gravitational collapse, scale echoing. 

The existence of this intermediate attractor gives rise to scaling relationships for near critical 
solutions. For any one parameter family of data in the limit of small black holes one finds a mass 
scaling relation: 

ln{MBH) lln\p - + f{ln\p - p^\), (1.2) 

where Mbh is the black hole mass on formation, 7 is a real parameter referred to as the critical 
exponent, f{ln\p — p.t\) is a periodic function, p is any parameter in the initial data and is the 
threshold value of that parameter that separates the data into two categories: supercritical (leading 
to black hole formation) and subcritical (leading to matter dispersion). Two separate analytical 
studies ([3], [1]) have shown that DSS results in the period of f{ln\p — p*\) being 

T = A/27. (1.3) 

Note that f{ln\p — p.t\) is often referred to as a "small oscillation" or "wiggle", but its explicit 
functional form must be determined numerically [3]. 

Garfinkle and Duncan have shown that a similar scaling relation exists for the maximum cur- 
vature at the origin for subcritical solutions [5] . For any one parameter family of data in the barely 
subcritical limit one finds: 

In (MAXt {i?(r = 0)}) ~ -2-i\Ti\p - p^\ + g{ln\p - p^\), (1.4) 
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where i?(r = 0) is the Ricci scalar at the origin and g{ln\p — p^\) is a periodic component, again 
with period T. This relation is similar to that of mass scaling except for the factor in front of 
In |p — I being —27 rather than 7. This derives from mass having units of length and curvature 
having units of length"-^ (with c = h = 1). 

The third aspect of critical collapse we note here is universality. DSS and critical scaling have 
been observed in the gravitational collapse of many different types of fields. The constants 7 and 
A are determined by the properties of the critical solution and are therefore universal: they vary 
for different kinds of matter but are independent of the form of initial data and the parameter p in 
that initial data that is varied. Numerical and semi-analytical calculations for the massless scalar 
field have given 7 ~ 0.37 and A ~ 3.4 ([3], [1], j6], [7], [5]). In all previous calculations, done in 
Schwarzschild or double null coordinates, the functions / and g were well approximated by a small 
amplitude sine function with period T = A/27 ~ 4.6. 

In this paper we study numerically black hole formation from the spherically symmetric collapse 
of a massless scalar field in Painleve-GuUstrand coordinates ([9], [l0])[^ PG coordinates were used 
by Wilczek and collaborators [Tl] to study Hawking radiation and more recently were proposed 
by Husain and Winkler [12] as a reasonable starting point to investigate the quantum dynamics 
of black holes. The key distinguishing feature of PG coordinates is that they are regular across 
apparent horizons. This allows us to evolve initial data until the onset of singularity formation 
when our chosen boundary conditions (zero mass density at the origin) are no longer valid. 

We first of all take advantage of the PG coordinate system to run the code beyond initial horizon 
formation to study in detail the formation of apparent horizons and the evolution of the resulting 
trapping surfaces. We then study the critical behaviour of the system and find our data to be in 
good agreement with previous results for the critical exponent and echoing period. Interestingly, we 
find the maximum curvature scaling relation to show the usual small oscillations but discover a new 
form for the oscillations in the mass scaling relation. The form of these mass scaling oscillations is 
found to be universal in the sense that it is the same for all families of initial data (that we tested) 
specified in PG coordinates. They are however quite different from those observed in Schwarzschild 
or double null coordinates. The reason for this difference is as yet unclear, but it is almost certainly 
related to the fact that the location of the initial apparent horizon can depend on the choice of 
spacetime slicings. 

While these results are interesting in themselves, our analysis is also important because it paves 
the way for studying the effects of quantum corrections on both the critical behaviour and the 
dynamics of trapping surfaces [13] . 

The paper is organized as follows: in Section [2] we lay out the foundation of our mathematics 
and develop the equations of motion used in our computer code. Section [3] describes our code 
giving relevant details of the numerical techniques and boundary conditions. Following this, we 
present the results of our simulations, discussing the scaling relationships and global properties of 
the spacetime (including trapping surfaces). We conclude with some remarks on these findings and 
consider future directions for this work. 

2 Equations of Motion 

The dynamical system that we wish to study is that of a collapsing spherically symmetric, massless 
scalar field in four spacetime dimensions. It turns out to be convenient to express these equations 
in a somewhat non-traditional parameterization, using the formalism of 2d dilaton gravity (see 

^ These coordinates are sometimes referred to as "fiat slice" coordinates because in vacuum tlie surfaces of constant 
time are flat. 
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for a general review of dilaton gravity and [15] for applications to black holes). In addition 
to providing a relatively simple form for the evolution equations in 4 dimensions, the formalism 
is generally applicable to a variety of spherically symmetric theories in four or more dimensions 
|8]. We will describe the formalism in some detail because the gauge fixing is a crucial part of 
our analysis. We begin by stating the classical action and defining a metric. We then perform a 
canonical transformation and fix the gauge to yield dynamical equations for the scalar field in PG 
coordinates. Throughout this paper we use a prime to represent dr and an overdot for dt- 
The action for a massless scalar field ip coupled to a dilaton in 2d is 



(2.1) 



where G is Newton's consant in 2d, g is the metric determinant, R{g) is the Ricci scalar and /i(0) 
is the dilaton coupling. 

The relationship of this action to spherically symmetric gravity in D = n + 2 spacetime dimen- 
sions is given as follows. The physical D metric is: 



dsl^^, = —^+r\c^)d^ll (2.2) 



where dO^ is the line element of the unit n-sphere and 



m := / d~W{^). (2.3) 







With the further substitutions 



^ .^5^, (.4) 



(n - l)z^W/"' 



(n-1) V/ 



(2.5) 



n \l . 

ijj = (2.^ 



the dilaton action (2.1 ) is (up to boundary terms) precisely equal to that of a spherically symmetric 



massless scalar field ifj minimally coupled to Einstein gravity in D spacetime dimensions: 

I(^) = j - \ I d^x^^)m\ (2.9) 

In the context of spherically symmetric gravity, performing the dimensional reduction in the action 
before varying yields the same dynamics as imposing spherical symmetric in the field equations 
themselves. 

In order to make our coordinate choice explicit and rigorous, we now summarize the Hamiltonian 
analysis for the theory. We first define the metric in modified ADM form: 

ds^ ■= e^P [-a^dt^ + {dr + Ndtf] (2.10) 
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where p, a and N are arbitrary functions of the spacetime coordinates. 

The authors of jl6| derive a partially gauge fixed Hamiltonian. For completeness, we review 
this procedure before applying a second gauge fixing condition that yields the completely gauge 
fixed metric in PG form. 

The first gauge choice is 

m = ^', (2.11) 



which yields the consistency condition on the lapse and shift functions: 

aGUp = N(t>'. 

The equations are put in a more transparent form using the canonical transformation: 



X :-- 
P :-- 



(2.12) 



(2.13) 
(2.14) 



where Hp is the momentum conjugate to p. P is then conjugate to X with their Poisson bracket 
being 

{X[x),P{y)] = G5{x,y). 



(2.15) 



The resulting Hamiltonian is 

i/(x,p,v,n^) 



dr a{ --^M' + Qm+1 



m' 



m 



where 



M 

Qm 



/-(tI-)' 



I {<!>'? I m 

2G\ X"^ P 



(2.16) 



(2.17) 
(2.18) 



M is called the mass function because it approaches a constant at spatial infinity where it is equal 
to the ADM mass of the solution. Qjn is the energy density of the scalar field. 
We now completely fix the gauge by choosing 



(2.19) 



This condition, along with Eqs.(2.17) and (2.11 ), implies that P = 2GM./1. These gauge conditions 



produce the non-static generalization of PG coordinates, as can be seen by using (2.12) and the 
gauge conditions in the case of a vacuum to reduce the line element to PG form |17] : 

/ I \ 2" 



-dr + \ dr + 



I2GMI 

~1W 



dt 



(2.20) 



More importantly, the spatial slices will be seen to be regular across apparent horizons that form 
during the evolution. 
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We are now able to write the equations of motion for the scalar field in fully reduced form: 



'lyj2GM/W 



+ 



n„ 



where a and ^A are the solutions to 



M' = gM + ^V^'n, 



V2GMIM 



I2GMI 



a 



(2.21) 



(2.22) 



(2.23) 



(2.24) 



The right hand side of (2.23) defines the instantaneous mass density of the configuration. It has 



the expected contribution from the energy density of the scalar and an additional coupling term 
that corresponds to the contribution from its self gravity. 

The above equations need to be supplemented by boundary conditions for P and a. Without 
loss of generality we choose cj = latr = 0. A change in this value corresponds to a trivial rescaling 
of the time coordinate. We also fix 7W = at r = 0, which guarantees that the metric is flat in 
the neighbourhood of the origin. In fact, a non-zero value of M at the origin signals the formation 
of a singularity, so that this boundary condition would have to be relaxed/modified in order to 
integrate the equations past singularity formation. 



3 Computational Methods 



Since we work in four dimensions (n=2), we have j{cl)) = 4' = r^/(4/^) and h{(j)) = Acj). For 
convenience we choose 2G = 1, which from Eq.( |2.4[ ) corresponds to ch oosing the characteristic 
length scale I to be Ipi = ^/G^. Note that the first gauge choice (2.11) implies that our spatial 
coordinate is the radius r of invariant two spheres. 

The first step in the iteration process is to specify initial ip and 11^ configurations. We work 
with two forms: 



Ar exp 



^tanh 



B 



r — vq 
B 



(3.1) 



(3.2) 



where A, B and ro are the parameters which may be varied to study critical collapse. We shall 
refer to these two forms of initial data as gaussian and tanh respectively. Note that since the mass 
density depends on ip' , the tanh form has one mass peak while the gaussian data has two. For both 
cases we define an initial standing wave by choosing n^(r, 0) = 0. 

The most interesting behaviour occurs as matter approaches the origin on ever finer scales. In 
order to observe this small r behavoiur while keeping a large enough lattice to contain all of the 
mass within the system (numerical errors occur when mass leaves the grid) , we refine the r-spacing 
Ar(r) near the origin so that it is several orders of magnitude smaller than the spacing near the 
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outer endpoint. We use an adaptive time step At{t) refinement according to the minimum found 
across the spacial slice using the condition 

At(t) = MIN, |^Ar(r)| (3.3) 

where ^ is the inverse of the local speed of an ingoing null geodesic. This provides stability by 
preventing information from moving over too many r-points in a single time step. Having At 
proportional to Ar greatly increases computational times for finer r-spacing, placing a practical 
lower limit on the size of black holes that we are able to create. However, this is not a significant 
concern as we are able to work on small enough scales to produce critical mass scaling. 

Both time and space integrations are performed using Runge-Kutta methods accurate to fourth 
order in the local grid spacing. For derivatives, we experimented with finite differences and cubic 
splines finding the same qualitative behaviour using both, which gives some indication of reliability. 
In the end we chose to use finite differences exclusively since they require less computational time. 



4 Results 



4.1 Black hole evolution 

Since our coordinates are regular across apparent horizons, we are able to evolve well beyond initial 
black hole formation up to the point when a singularity begins to form and our boundary conditions 
at the origin are no longer valid. Spacetime diagrams of black hole formation and evolution are 
given in figures 1(a) [ and 1(b) for tanh and gaussian initial data respectively. Using this same data, 
we have also constructed animations of the mass density profile throughout this process, available 
for viewing at http:/ / theoryxS. uwinnipeg. ca/ users/ jziprick/. 




space space 

(a) tanh initial data (b) gaussian initial data 



Figure 1 : Spacetime diagrams of the collapsing massless scalar field using radius for the space axis 
and PG time for the time axis. The solid black lines are null geodesies and the dashed red lines 
indicate the trapping surface boundary. 

We note that apparent horizons form in pairs within constant PG time slices. This is to be 
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expected from looking at the apparent horizon condition: 

\Mh? = 0, (4.1) 
^rh-2G^^^M = 0, (4.2) 

which is initially positive across the r-axis. Minima begin to form as mass collapses inward. Black 
hole formation is indicated by the appearance of a single horizon where the function first dips to 
zero. This horizon splits into two as the minimum drops below zero; the outer horizon continues to 
grow as mass passes through it while the inner horizon moves toward the origin. Depending upon 
the initial mass distribution, subsequent horizon pairs may form as demonstrated in figure |l(b) 
which depicts the evolution of gaussian initial data. The horizon paths define a two dimensional 
trapping surface within which all null geodesies point towards the origin. 



4.2 Critical behaviour 

In our study of the critical behaviour, we varied the parameters B, and in the gaussian form of 
initial data, and A in the tanh form. The critical values were found to an accuracy of dp/p ~ 10^^^ 
using a binary search. We tried a variety of grid resolutions with r-spacings near the origin ranging 
from 10" 5 to 10^^, which in turn affects the t-spacings. The accuracy was roughly equal in all 
cases where black hole size was not significantly affected by the resolution (when r/j >> Ar(r/i)). 
However, finer grids allowed for data to be found for lower values of In^p — p^,) thus showing critical 
scaling over a wider range of a given parameter. 

To demonstrate scale echoing, we present a plot of the scalar field at r = as a function of 
r = /n(t* —t) for the nearest solution to criticality that was numerically achieved using an r-spacing 
near the origin of 10"^. t=K is the time of black hole formation in the critical solution given by the 
value approached asymptotically in the small mass limit. Since we choose a = 1 as our boundary 
condition at r = 0, t is equivalent to the proper time at the origin. ip(r = 0) is seen to be periodic 
in T with a period of A = 3.43 it 0.06 found by averaging over the two oscillations between the 
minima at r ~ —2.5 and —9.3, with an error determined by the uncertainty in t*. This result is in 
good agreement with the accepted value of A ~ 3.44. 

Table [T] lists the critical exponent and period found for various lattice resolutions and families 
of initial data using both the mass and maximum curvature scaling relations. Our results in all 
cases are seen to agree with the accepted values of 7 ~ 0.374 and T ~ 4.6. 



type of 
scaling 


form of 
initial data 


parameter 
varied 


resolution 
near origin 


7 


T 


curvature 


gaussian 


A 


10-5 


0.375 ± 0.005 


4.61 ±0.05 


mass 


gaussian 


A 


10-5 


0.375 ± 0.003 


4.6 ±0.1 


mass 


gaussian 


A 


10-4 


0.374 ± 0.003 


4.4 ±0.2 


mass 


gaussian 


B 


10-4 


0.376 ±0.006 


4.4 ±0.2 


mass 


gaussian 


ro 


10-4 


0.377 ±0.005 


4.4 ±0.2 


mass 


tanh 


A 


10-4 


0.371 ± 0.004 


4.5 ±0.1 



Table 1: The critical exponent and period of the mass scaling curve for various initial data and 
mesh resolutions. 



A plot of the maximum subcritical curvature scaling behaviour is given in figure 3(a) 



7 was 

determined using a least squares linear fit over the points where the solution is seen to follow the 
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Figure 2: A plot of the scalar field at the origin versus r illustrating scale echoing in the critical 
solution. The r-spacing in generating this data was 10~^ near the origin. 



intermediate attractor. The error quoted indicates how this value changes depending on the range 
of points chosen. The period and associated error (given in the first row of table [T]) were determined 
by fitting a sine function to the residual left after subtracting the linear fit from this data. 



Plots of the supercritical mass scaling behaviour are presented in figures 3(b) -3(d) On each of 



these we have drawn a line that osculates the peaks of the curves. The values of 7 were determined 
by using a least squares linear fit through the data points that are very near this line (within ~ 5 % 
of Mb//)[^ The error quoted is twice the standard deviation in the calculated slope. The periods 
are found by averaging the distance between the cusp-like valleys, with the error representing the 
uncertainty in locating the cusps. The data used in these calcultions consisted of the two periods 
corresponding to the lowest range of the parameter that produced black holes such that the size was 
not significantly affected by lattice resolution (for example, the periods between /n|A — ~ —9 
and —18 in figure 3(c) were used to calculate the values given in the second row of table 

Previous studies of critical collapse (done in either Schwarzschild or null coordinates) have shown 
both scaling relationships to be a slight oscillation about a straight line that is well approximated 
by a sine function. While this is confirmed by our results for the subcritical maximum curvature 
scaling, it is clearly not the case for our supercritical mass scaling results where we find a periodic 
form with relatively large amplitude and a distinctly different functional form. Considering that 
our results agree with the known parameters of Choptuik scaling for varied initial data and mesh 
resolutions, and that our new-found periodic form is universal within the class of coordinate systems 
that we are using, it is unlikely that this periodic form is a numerical artifact. 

According to Gundlach [3_ the function / that determines the small wiggle is universal. However, 
the location and time of formation of the initial apparent horizon can depend on the choice of 
spacetime slicing. Ours is the first calculation in PG coordinates. These are non-static coordinates 
for which data are specified on spacelike surfaces that are regular across future horizons. These 
unique features of PG coordinates are undoubtedly at the root of the new periodic form in the 
mass scaling relationship. It would be of great interest to continue the near critical evolution past 
initial horizon formation in order to study explicitly the scaling relation of the final horizon radius. 



^Note that our results are insensitive to this choice at the given level of precision. 
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which is independent of spacetime slicing. This is currently under investigation. 

5 Conclusion 

We have successfully demonstrated Choptuik scaling for a massless scalar field in PG coordinates 
and revealed an interesting large amplitude component to the oscillations in the mass scaling 
relation. Our gauge choice allows integration of the dynamical equations beyond horizon formation 
until the onset of singularity formation. This allowed us to study in detail the evolution of black 
holes showing explicitly that apparent horizons form in pairs with the number of pairs depending 
on the complexity of initial data. 

One possible direction for work with this system would be to take the number of dimensions as 
a parameter and compare the results with the dimensional analyses done in [8] and p8l . 

Particularly exciting is the prospect of incorporating quantum corrections into the gravitational 
potential^ Loop Quantum Gravity suggests that the underlying discreteness at the quantum level 
will give rise in the semiclassical limit to an effective short range repulsive component to the 
gravitational potential which in turn can lead to singularity avoidance |19| 1^ . Husain has 
done an investigation of the effect of such a correction on Choptuik scaling using double null 
coordinates which did not allow the code to be extended past horizon formation. Since we are 
able to investigate the dynamics across the entire spacetime we can investigate directly the effects 
of short range modifications to the gravitational potential near the onset of singularity formation. 
We have begun such an investigation |13j and preliminary results indicate dynamical singularity 
avoidance that allows us to calculate numerically the non-singular spacetime showing black hole 
formation and evaporation. These results will be presented in a future publication. 
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ln(MAX,{R(r=0)}) = -0.7528 ln|A-A.| + 2.637 




ln|A-A.| 



(a) The maximum curvature scaling relation for a 
variation in the paramter A in gaussian initial data 
with Ar = 10~^ near the origin. 




ln|A-A.| 



(b) The mass scaling relation for a variation in the 
paramter A in gaussian initial data with Ar = 10~* 
near the origin. 




n|A-A.| 



(c) The mass scaling relation for a variation in the 
paramter A in gaussian initial data with Ar = 10'^ 
near the origin. 




ln|A-A.| 



(d) The mass scaling relation for a variation in the 
paramter A in tanh initial data with Ar = 10~^ near 
the origin. 



Figure 3: Plots of the critical scaling behaviour. 
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